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GLOBAL DISCRETE ARTIFICIAL BOUNDARY CONDITIONS FOR 
TIME-DEPENDENT WAVE PROPAGATION* 

V. S. RYABEN’KIlt, S. V. TSYNKOVS* AND V. I. TURCHANINOVN 

Abstract. We construct global artificial boundary conditions (ABCs) for the numerical simulation of 
wave processes on unbounded domains using a special non-deteriorating algorithm that has been developed 
previously for the long-term computation of wave-radiation solutions. The ABCs are obtained directly for the 
discrete formulation of the problem; in so doing, neither a rational approximation of “non-reflecting kernels,’ 1 
nor discretization of the continuous boundary conditions is required. The extent of temporal nonlocality of 
the new ABCs appears fixed and limited; in addition, the ABCs can handle artificial boundaries of irregular 
shape on regular grids with no fitting/adaptation needed and no accuracy loss induced. 

The non-deteriorating algorithm, which is the core of the new ABCs, is inherently three-dimensional, it 
guarantees temporally uniform grid convergence of the solution driven by a continuously operating source 
on arbitrarily long time intervals, and provides unimprovable linear computational complexity with respect 
to the grid dimension. The algorithm is based on the presence of lacunae, i.e., aft fronts of the waves, 
in wave-type solutions in odd-dirnension spaces. It can, in fact, be built as a modification on top of any 
consistent and stable finite-difference scheme, making its grid convergence uniform in time and at the same 
time keeping the rate of convergence the same as that of the non-modified scheme. 

In the paper, we delineate the construction of the global lacunae- based ABCs in the framework of a 
discretized w r ave equation. The ABCs are obtained for the most general formulation of the problem that 
involves radiation of waves by moving sources (e.g., radiation of acoustic w r aves by a maneuvering aircraft). 
We also present systematic numerical results that corroborate the theoretical design properties of the ABCs’ 
algorithm. 

Key words, artificial boundary conditions, wave propagation, lacunae, non-deteriorating method 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Numerical simulation of wave phenomena on unbounded domains (e.g., radiation 
and/or scattering of acoustic or electromagnetic w r aves), often encounters the following well-recognized dif- 
ficulty. As no computer can either handle infinite arrays of data or perform infinite numbers of arithmetics 
operations, the discrete approximation of the problem has to be made finite (i.e., finite-dimensional). Con- 
sequently, the original infinite domain has to be truncated and special artificial boundary conditions (ABCs) 
have to be developed as a closure for the resulting finite formulation. 

The literature on the subject of ABCs is very extensive, see, e.g., review" papers by Givoli [1] and 
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Tsynkov [2], as well as another recent review by Hagstrom [3], which is geared more specifically toward wave 
propagation problems. In the current study, we are going to focus on genuinely unsteady (as opposed to 
time-harmonic) wave phenomena to be computed in the time domain. For this type of problems, most of 
the ABCs’ research to date has been done in the framework of simple approximate local methods based, 
e.g., on quasi-one-dimensional characteristics’ decomposition. These methods often appear insufficiently 
accurate. Some of the more accurate methods that have recently gained attention are based on the so-called 
perfectly matched layers (PML), see the original publications [4-7] and reviews [8,9]. Unfortunately, as 
shown in [10, 11], these methods may give rise to instabilities of a particular kind. The latter typically 
manifest themselves when integrating over long time intervals and thus exacerbate even further the well- 
known problem of accumulation of error in long-term numerical simulations. 

Among other existing unsteady ABCs’ approaches, only very few methodologies can guarantee the accu- 
racy that theoretically would not hamper that of the interior approximation. All of these methodologies are 
non-local, see, e.g., [12 18], which is characteristic of highly-accurate (ideally, exact) ABCs. The techniques 
of this group typically involve two “approximating” steps, which are undertaken consecutively when build- 
ing the ABCs. The first step is the replacement of the fully non-local in space- time true exact boundary 
conditions, which are most often written using pseudo-differential operators, with approximate boundary 
conditions (still at the continuous level) that would provide for only a limited extent of non-locality. More 
precisely, this step aims at limiting the temporal nonlocality of the ABCs, which may be prohibitively expen- 
sive in computations. This can be achieved, e.g., by employing a rational approximation of the corresponding 
symbol (kernel). 1 The first step is then followed by the second one, which is the discretization of the result- 
ing localized continuous boundary conditions. We note that unless given a special thorough attention, the 
discretization step may lead to accuracy loss and again cause instability (this pertains to purely local ABCs 
as well). We also note that all of these techniques are restricted geometrically to computational domains of 
a simple regular shape, e.g., those with spherical or planar boundaries. 

Our recent w T ork [22,23] indicates that the issue of time-dependent ABCs may be closely related to 
another problem that has been mentioned before — the accumulation of numerical error during long runs. 
This problem has been recognized as an outstanding question in computational PDEs for many years, since 
the first systematic convergence studies for discrete approximations have been conducted in the fifties. From 
the standpoint of practical computing, deterioration of the solution over long time intervals can be attributed, 
e.g., to the mechanism of either numerical dissipation or dispersion or both. Theoretically, this phenomenon 
is conveniently termed as non-uniformity of the grid convergence in time, and all conventional discrete 
approximations that are currently in use are known to suffer from this deficiency. 

As our w T ork [22, 23] suggests, the key to resolving the question of long-term error accumulation may 
be in precisely following the physical nature of the original problem when building a numerical algorithm. 
Namely, it is known that waves in three dimensions have aft (or trailing) fronts. This is a manifestation of 
the so-called Huygens’ principle or more generally, the presence of lacunae in w T ave-type solutions in odd- 
dimension spaces. Using this property, we have been able to develop a modification to any consistent and 
stable finite-difference scheme that approximates a Cauchy problem for the wave equation so as to make its 
grid convergence uniform in time on arbitrarily long intervals. The uniform temporal convergence of the 
algorithm has been proven theoretically along wdth its optimal computational complexity (i.e., linear with 

! In fact, a wide variety of purely local ABCs (i.e., local in both space and time) can be obtained via rational approximation 
of symbols as well; this approach has been known for two decades, see [19-21]. The general trend is that the more of the 
nonlocal nature of exact ABCs is compromised, the less accuracy one can expect from the resulting approximate methodology. 
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respect to the grid dimension). The rate of temporally uniform grid convergence, see [22,23], remains the 
same as that of the original non-modified scheme. These results apply to the general case of moving sources 
of waves that operate continuously in time. As an example, we show in [23] that the linearized flow around 
a maneuvering aircraft can be interpreted in this framework. 

At least as important, the procedure of [22,23] allows one to replace the original infinite domain of the 
initial-value problem by a finite computational domain that would facilitate the construction of a finite- 
dimensional discretization. As will be seen from the discussion in the current paper, the latter replacement 
leads to obtaining highly accurate non-local unsteady ABCs for combined problems that may include complex 
phenomena on a bounded interior domain but reduce to the homogeneous wave equation in the far field. 
Similarly to the case analyzed in [22,23] the interior domain may be moving. The ABCs are built directly 
for the specific interior approximation used and in this sense can be considered its most natural extension. 
We emphasize that they involve neither of the two common approximating steps (rational approximation 
of the symbol followed by discretization) that have been mentioned before in connection to some existing 
ABCs’ methodologies. 

Unlike in all other methods available in the literature, the extent of temporal nonlocality of the unsteady 
ABCs based on the technique [22,23] is bounded naturally for all times because of the explicit use of lacunae. 
This bound is not a consequence of any approximation. Moreover, as these ABCs are obtained directly for 
the specific finite-difference scheme, the issue of discretization of the boundary conditions, which has been 
shown to cause problems before, simply does not arise in this framework. Besides, the new ABCs possess 
full geometric universality, i.e., can handle any shape of the external artificial boundary with equal ease on 
a regular Cartesian grid with no fitting/adaptation required and no accuracy loss caused. 

The rest of the paper is organized as follows. In Section 2, we provide for a brief outline of the phe- 
nomenon of lacunae in wave radiation solutions, and show how one can use those to obtain a non-deteriorating 
algorithm for the long-term numerical integration of the corresponding problems. In addition to the theoret- 
ical justification, we also include in this section several computational demonstrations of the properties of the 
aforementioned algorithm. In Section 3, we describe in detail the construction of the global finite-difference 
lacunae-based ABCs, and briefly comment on how the proposed construction fits into the general framework 
of discrete time-dependent boundary conditions developed by Ryaben’kii in [24]. Finally, Section 4 contains 
an extensive set of numerical experiments with the new ABCs for the wave equation. The experiments 
are conducted for finite- difference schemes of different orders of accuracy, different laws of motion for the 
waves 1 sources (uniform, as well as non-uniform), and different interior models that require closure by the 
homogeneous wave equation in the far field. These experiments corroborate the theoretical design properties 
of the ABCs 1 algorithm. 
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2. Lacunae and Non-Deteriorating Numerical Integration. 


2.1. Lacunae of the Wave Equation. We consider a Cauchy (initial- value) problem for the three- 
dimensional wave equation, x = (xi , £2,2:3): 


gv _ -2 fgV gv gv \ 

\ 3x 2 dx. 2 <9x| / 


/(*,<), i>0, 


(2.1) 


V 5 


<=0 



= 0 . 


( 2 . 2 ) 


(The limitation of having homogeneous initial conditions (2.2) can be alleviated, see [22,23].) The problem 
(2.1), (2.2) is solved on the domain S(t) C M 3 , which has finite diameter d for all times t > 0; other than that 
the domain S(t) may travel in space according to an arbitrary law of motion except that its maximum speed 
k is required to be “subsonic:” k < c. The solution <p{x, t ) is driven by the continuously operating source 
/(x,£), /(x, 0) = 0, and we require that Vt > 0: supp f{x y t) C S(t). In other words, we study the radiation 
of waves by a source, which is compactly supported in space for all times. The solution is of interest for us 
also on a compact domain, which we call 5(t); it fully contains the source and follows its motion if there is 
motion. This is a simplified model for many interesting physical phenomena that are more complex in their 
nature. As we shall see, this model is very useful when constructing the ABCs for a variety of problems. 

For every (a?, f), the solution ip = p(x, t) of problem (2.1), (2.2) is given by the Kirchhoff integral: 


ip(x, t) 



e<ct 





(2.3) 


where £ - (6,6.6), Q - \x - £| = \A z i ~ 6) 2 + (*2 - 6) 2 + (^3 - 6) 2 , and = d^d^d^- ^ we 
assume for a moment that the right-hand side (RHS) f{x,t) is compactly supported in space-time on the 
domain Q C E 3 x [0, -foo), then formula (2.3) immediately implies that 


ip{x, t) = 0 for (x, t) € f| {(*,t)|l*-*l <«#-*), t >e). (2.4) 

(£J)eQ 

The region of space-time defined by formula (2.4) is called lacuna of the solution <p = ip(x,t). This region 
is obviously obtained as the intersection of characteristic cones of equation (2.1) once the vertex of the cone 
sweeps the support of the RHS: supp / C Q. From the standpoint of physics, the lacuna represents that part 
of space-time where the waves generated by the sources /, supp / C Q , have already passed and the solution 
has become zero again. (Sometimes, the name “secondary lacuna” is used to distinguish it from the primary 
lacuna, which is the area where the waves have not reached yet.) The phenomenon of lacunae is inherently 
three-dimensional. The interior surface of the lacuna represents the trajectory of aft (trailing) fronts of the 
waves. The presence of aft fronts in odd-dimension spaces is known as the Huygens 1 principle, as opposed 
to the so-called wave diffusion which takes place in even-dimension spaces. 

2.2. Computation with a Compactly Supported Source. Assume now that the moment of in- 
ception of the source /(x,f) is to (in particular it may be t 0 = 0); at this moment the domain S(to) of the 
RHS /(x, t) occupies a position in space which is schematically represented by the interval [A\, A 2 ] of size d 
on Figure 2.1. 2 Assume also that by the time t\ > to this source ceases to operate, which makes the RHS of 

throughout this section we will be using schematic one-dimensional illustrations always keeping in mind, however, that 
the actual problem is three-dimensional. 
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Fig. 2.1. ID schematic representation for the fronts of waves generated by a compactly supported source. 

equation (2.1) compactly supported in both space and time: supp / C Q — {(x,t)|x € S(t ), to < t < ti] . 
Clearly, by the time t \ the domain S(ti ) can be displaced from its initial location no further than the distance 
k(t\ — to) is each direction, which is schematically represented on Figure 2.1 by the boundaries of the interval 
[Bi y B 2 ] of size d + 2k(ti — t 0 ). Starting from t — t\ no new waves will be generated, and those generated 
prior to t\ will continue traveling in space and thus will eventually leave the domain S(t) completely, because 
their speed of propagation c is higher than that of the domain, k. The moment t 2 when this happens, i.e., 
when the solution (p(x,t ) again becomes zero on 5(t), is easy to calculate, see Figure 2.1. By this moment, 
the domain S(t 2 ) can travel no further than the interval [C\,C 2 ] of size d + 2k(t 2 — fo)> and we need to 
assume that the aft fronts will also be exactly at the boundaries of this interval at t — t 2 , which immediately 
yields 

= 0, for x e 5(t), t > t 2 = t 0 + + k) ^ ^2.5) 

c — k 

Estimate (2.5) is fundamental. It essentially says that once w r e need to calculate the solution ip (x,t) on S(t) 
and the sources are compactly supported in space-time: supp / C Q { (x,<)| x E 5(£), to < t < then we 
may stop the calculation at t — t 2 because afterwards the solution on S(£) will be zero anyway. This means, 
in particular, that if the solution is calculated using a discrete method, e.g., a finite-difference scheme, then 
no new error will be accumulated after t — t 2 . The constants in both consistency and stability estimates 
of the scheme (see [22,23] and below r for detail) will now depend on the time interval T mt = 
rather than final time Tfi na i, which for the case of a compactly supported RHS simply becomes immaterial. 

Besides, once we stop the calculation at t 2 ~ t 0 4- T[ nt , we realize that during the time interval Xi nt that 
has passed since the beginning t ~ to, no waves could have traveled in space further than the boundaries 
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of the interval [Di,D 2 ] of size d + 2cT mU see Figure 2.1. Beyond this region the solution is zero because 
this is the area of the primary lacuna. Therefore, even though the original problem was formulated on an 
infinite domain, we can, in fact, calculate the solution on a finite domain [D \ , D 2 ] of size d + 2cT\ nt with zero 
external boundary conditions (of the Dirichlet type). 

The transition from the infinite domain to a finite one does not, obviously, come “at no charge.” One 
can rather say that it comes at the expense of having the computational domain [Di,D 2 ] larger than the 
actual domain of interest S(t). However, the size of the “redundant” portion of [Di,D 2 \ can be further 
reduced by observing that all we have to do is make sure that by t = t 2l i.e., by the moment the last waves 
generated by /(x,£), supp / C Q ) leave S(t), no new waves can enter S(0* This can be guaranteed either 
as indicated previously, by placing the outer boundary sufficiently far away so that no waves from /(x,0> 
supp / C Q, can even reach it until t = t 2 = £o + Xj nt or alternatively, by placing it closer and thus allowing 
for reflections, but still not too close so that no reflected waves can come back, i.e., reach S(t), by t = t 2 . 
The size of the new, smaller, computational domain [Ex,E 2 ] with reflecting outer boundary, see Figure 2.1, 
can be estimated easily. The minimum size Z, see Figure 2.1, is found by requiring that the reflected waves, 
wffiich travel with the same speed c but in the opposite direction, reach the boundary of [C\,C 2 ], he., the 
utmost possible location of S(t 2 ), by the exact same moment of time t — t 2 when the aft fronts leave S(t). 
This immediately yields 

Z = d+(* + c)r int . (2.6) 

By comparing the value of Z from (2.6) with the size of jXb, £> 2 ]? which is d + 2 cXi nt , we conclude that the 
extra size of the computational domain beyond d can be reduced by up to a factor of 2 (when k = 0) in each 
coordinate direction. 

We also note that in fact any well-posed boundary condition can be specified at the reflecting outer 
boundary of [E\,E 2 ]. The reason is that this boundary is intentionally positioned so that the reflections are 
not going to have any effect on the solution inside S(t) anyway. A particularly convenient way to treat the 
boundary of [E\ , E 2 ] will be to set the periodic boundary conditions there. In so doing the three-dimensional 
rectangular domain becomes a three-dimensional toroidal surface (the opposite faces of the rectangle are 
identified with one another) and we only have to keep in mind that the reflected waves will now need to 
be interpreted as those that leave the domain on one side and enter it from the opposite side. This new' 
interpretation obviously brings no change to the foregoing considerations that led to the size estimate (2.6). 
However, for the case of a continuously operating traveling source that we analyze be!ow% periodicity implies 
that the motion of the source can also be formally considered on the toroidal surface, wdiich makes the 
computational setup much simpler. 

2.3. Computation with a Continuously Operating Source. Both foregoing observations — finite 
time interval T mt and finite spatial domain [E\^E 2 ] needed for calculating the solution driven by the sources 
supp / C Q = {(x,0| x E S(0> t Q < t < t\ }, on the domain S(t) — are crucial for the original case 
of a continuously operating source supp / C {(x,t) | x E S(t ), t > 0}. In this case we first take a 

parameter T > 0 and introduce a smooth even compactly supported function ©(/), t E K, of a “hat” type: 

0(t) = 0 , \t\ > T , 

m = ©(-*) , 

0(0 = 1 , t E [~oT , crT] , 0 < a < 1 , 

e (H £r+< ) = 1 - e ( 1 r £T - 1 ) ■ 
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which obviously generates a partition of unity: 


l = ^e(*-(l+<r)Tj) , < > 0 , 

3=0 

with the overlap size (1 — a)T. Then, we represent the right-hand side of equation (2.1) in the form 

oo oo oo 

/(*, t) = /(*, t ) £ 0(f - (1 + a)Tj) = Y, 0(t - (1 + o)Tj)f[x, <) = £ /#(*. t, r) • (2-7) 

>=o i— o j=o 


Clearly, for each fj(x,t,T) = 0(# — (1 + cr)Tj)f(x , £), j =0,1,..., we have 

supp/j(x,£,T) C {(x,£)|x € 5(<), (j - 1 + ^j)r < t < (j + 1 + vj)T} . (2.8) 


Due to the linear superposition, the overall solution <p(x,t) of equation (2.1) will be given b}' the sum of 
individual contributions from fj(x, t, T), j = 0, 1, . . . : 

oo 

¥>(*. *) = l ' T ) > ( 2 - 9 ) 

j=o 


where each contribution ipj(x,t>T) solves the following sub-problem: 

gVj _ r 2 f \ — f ix t T\ 
df- \ dx\ + dx\ + dx 2 3 ) ~ ’ ,T) 1 

_ Ml 


V>i 




dt 


( 2 . 10 ) 


= 0, J = 0, 1,2, . . . 


t=(i-l+<rj)T 


Notice that each ^ (x, t,T), j = 0, 1, . . . , can be calculated absolutely independently of the others and that 
the corresponding source term fj(x, t, T) is a function compactly supported in both space and time, see (2.8). 
Consequently, according to (2.5), if we interpret t 0 and t\ as (j — 1 + a j)T and ( j -b 1 + oj)T y respectively 
(see Figure 2.1), then we can conclude that every (pj(x y t,T ) of (2.9) needs to be calculated only during a 
finite interval of time T mt = . It is important to realize that this interval does not depend on the 

actual moment of time t. 

Moreover, even so the series (2.9) is formally infinite, it is easy to see that for any t > 0, x £ S(t ), it 
contains only a finite fixed number of non-zero terms. First of all, because of the causality, <^j(x,t,T) — 0 
for x £ S(t) if t < (j — 1 + crj)T. In other words, for a given moment of time t, the contribution of all those 
fj(x,t,T) that are active only at subsequent moments of time, is obviously zero. A somewhat less trivial 
observation is that because of the lacunae the contribution of the “sufficiently retarded” terms /j(x,t,T) to 
the overall solution at a given time level t will be zero as well. More precisely, <^(x,t,T) = 0 for x € S(t) 
if (j — 1 + aj)T < t — Tj nt . This follows immediately from (2.5) assuming that to = (j — 1 + aj)T and 
t\ = (j + 1 + crj)T . Consequently, instead of (2.9) we can write: 


if{x,t) = , x 6 S(t) , 


( 2 . 11 ) 


where p\ = -b l)J , pi = + l) j , and [ • ] stands for the integer part. The expressions for 

Pi and p 2 indicate that we will always have either p\ = p 2 - [ (i+£)t ] or P l ~ P 2 “* ] — T Therefore, 

the number of terms p = p 2 — (pi — 1) in the sum (2.11) will never exceed + 2. As T int does not 
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depend on t we conclude that neither does the foregoing upper bound for p . As such, the number of terms 
in the sum (2.11) can always be considered finite and fixed. Altogether we obtain that Vt > 0 the solution 
x E S(t), is composed of a finite non-increasing number p or additive terms, and each of the latter 
needs to be taken into account only during a finite non-increasing interval of time T \ nt . 

In the perspective of numerical computation, the latter consideration translates into temporal^ uniform 
grid convergence of the discrete algorithm. Indeed, assume that we are integrating equation (2.1) by means 
of a finite-difference scheme with the order of accuracy 0(h a ), where h is a general notion for the grid size 
and a > 0. Then, the discrete solution <p^(x,£) converges to the continuous solution as the grid 

size decreases: 

-¥>(*, i)|| <K h a , t e [0,T final ] , (2.12) 

here Tfi na i is the total integration time. Inequality (2.12) is a generic convergence estimate; it holds provided 
that the RHS /(x, t) of equation (2.1) is sufficiently smooth. A detailed discussion on the smoothness require- 
ments for f(x,t) can be found in [23], along with the specific consistency /stability /convergence estimates in 
the norms that would take into account a particular smoothness level. 3 

The constant K in inequality (2.12) does not depend on the grid. It is, however, knowm to depend on 
the actual RHS /(x,£), as well as the final time T^ na \: K = K(f, 7k n ai)- The dependency of K on Tfj na \ is 
typically a growth, and sometimes this growth may be rapid. This means that even so on any fixed interval 
[0 : Tfi na i] the scheme converges as h — » 0, to obtain the same level of accuracy on a larger [0, Tfi na i] one may 
need to take a finer overall grid ahead of time. Thus, the convergence appears temporally non-uniform. On 
the language of practical computing, this phenomenon can be interpreted as the accumulation of numerical 
error over long runs. This issue has been long acknowledged unresolved in the literature. 

The situation changes dramatically if, instead of the straightforward time-dependent integration, we 
first use the foregoing lacunae-based representation (2.11) of the solution <p(x,£). In so doing, for each 
j == pi,... ,j> 2 i we still integrate the corresponding sub-problem (2.10) using the same finite-difference 
scheme as before. However, convergence estimate for the scheme then becomes: 

x € S(t) , t e[(j -1+ aj)T } (j - 1 + (jj)T + Tint] . 

A very important circumstance is that unlike K in estimate (2.12), the constant K 3 in (2.13) for each j 
depends on T\ ni rather than Tf\ na j: Kj = Kj(fj 7 T ini ). Keeping in mind that each <p^ (x, t, T) can be 
computed independently of the others, and using linear superposition (formula (2.10)), we then easily obtain 
instead of (2.12): 

<p^(x,£) - <p(x,t) | < p • K • h a , x 6 S(£), t > 0 , (2.14) 

where K — K(f , T mt ). Note that p is fixed and does not increase with £, and K now depends on Tj nt rather 
than Tfinai? where T mt is also fixed and does not increase with t. Therefore, estimate (2.14) implies temporally 
uniform grid convergence of the discrete lacunae-based algorithm on arbitrarily long time intervals or in other 
words, for t > 0. A detailed formal proof of this result that, again, involves specific norms, can be found 
in [23]. Here w r e only need to add that for each inequality (2.13) to hold, the corresponding fj(x,t,T), 

3 Smoothness of the source terms will also be important when constructing the lacunae-based ABCs, see Sections 3 and 4. 


< Kj • h a 


(2.13) 
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j = p lt . . . ,p 2 , has to possess the same regularity as that required for the original scheme to converge. This 
explains why choosing the partition (2.7) smooth with overlaps was very important. 

From the standpoint of practical computing, temporally uniform grid convergence implies that the 
numerical error will not get accumulated beyond some predetermined bound for as long as the computation 
needs to be performed; and once the grid is refined the aforementioned bound will also drop in accordance with 
the specific rate 0(h Q ). This is clear because when the calculation is stopped for a given term 
after the interval T\ ni has elapsed, the error will not be accumulated any further, and the number of terms p 
that need to be taken into account is fixed and non-increasing. Thus, we have obtained a non-deteriorating 
numerical algorithm for integration of the wave equation over arbitrarily long times. Let us emphasize that 
it can be built as a modification of any consistent and stable finite-difference scheme, and that it preserves 
the original rate of convergence of the scheme while making the convergence uniform in time. 

Besides, let us assume, for example, that the original finite-difference scheme has linear computational 
complexity with respect to the grid dimension, which is typical for explicit schemes. Then, it is easy to see 
that the modified lacunae-based algorithm will also have linear computational complexity with respect to the 
grid dimension. Indeed, this immediately follows from the fact that each term t , T) is computed using 

the original scheme on a compact domain of size Z (see Figure 2.1) during a finite fixed interval of time 7j nt , 
and the number of terms p is, again, fixed and non-increasing. We should note that for the type of problems 
that we are studying linear complexity with respect to the grid is, in fact, optimal, i.e., unimprovable. 

2.4. Computation Using Continuous Time Marching. The following, and last, step in building 
the lacunae-based algorithm for long-term numerical integration of the wave equation is to realize that for 
implementing formula (2.11) we do not necessarily need to compute each term independently of 

the others. Instead, we can implement the algorithm in a way similar to the standard time-marching by 
means of a finite-difference scheme. For that, we will need to use the aforementioned periodic boundary 
conditions on the outer boundaries of the auxiliary domain [£i, £ 2 ], see end of Section 2.2 and Figure 2.1. 

The first key observation that we make here is that once the motion of the wave sources, as well as 
the propagation of waves themselves, are considered on a three-dimensional toroidal surface, rather than on 
the genuine R 3 , then for every portion of the RHS fj(x,t,T) it does not really matter where on the period 
this source is located, or where it starts its motion from, at to — (j — 1 + crj)T. It does not have to be 
exactly “in the middle” as shown on Figure 2.1, because all locations on the period (i.e., toroidal surface) 
are equivalent. All we have to worry about is that by the time £2 = (j — 1 + <rj)T 4- T \ n t the waves generated 
by fj(x,t,T ), see (2.8), will have left the domain S(t), and that no waves could have re-entered this domain 
during And this will be exactly the case because the size Z = d + (k + c)Tj nt , see (2.6), has been 

chosen sufficiently large to provide for that. Since we always assume (for simplicity) that the period Z is 
the same in all coordinate directions, then we only need to formally consider /j(x, t } T) instead of fj( x, t, T ) 
and accordingly, (pj(x,t, T) instead of T), where x = (xi>£ 2 >£ 3 )> and xi = xi — [^-] Z, i = 1,2,3. 

Next, we shall analyze formula (2.11) from a slightly different point of view. To begin with, we notice 
that on the initial stage of computation, i.e., when t is small, the lower summation limit p\ may turn out 
negative. Basically, it does not create any inconsistency and does not cause any problem because f(x y t) =0 
for t < 0 anyway. In fact, we can simply disregard all negative j ' s in the sum (2.11) for small t*s and 
initially consider the summation 0 (pj(x,t,T) instead of (2.11). “Initially” here means till the actual 
expression pi = -I- l)j becomes positive. It is easy to see that the computation on this initial 

stage is equivalent to the conventional time-marching of the wave equation (2.1) on the domain [Ei,Eo] of 
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size Z (see Figure 2.1) with periodic boundary conditions. Indeed, all we do here is simply take into account 
one component of the source fj after another. Due to linear superposition, this amounts to the continuous 
integration of the wave equation driven by / = fj fr° m £ = 0 till the actual time t. We also note that 
the duration of the initial stage is, obviously, T\ nt . And the period Z, see (2.6), has been chosen sufficiently 
large so that for the time interval of length 2] nt there will be no difference on the domain S(t) between the 
solution <p(x,t) computed in the periodic setting and the solution that one could have possibly computed 
with no periodization (see Section 2.2). 

As soon as the time interval Tint, = has elapsed since the inception moment t — 0, the 

computation enters its regular (as opposed to initial) stage. This regular stage, which can, in fact, be 
continued for as long as necessary, is characterized by the positive values of pi (the first positive value is 
obviously p\ = 1) and finite non-growing number p — p 2 — (pi — 1) of terms in the sum (2.11). 

On the regular stage of the algorithm, we continue marching equation (2.1) with periodic boundary 
conditions in space. Obviously, as the time t elapses both pi and in formula (2.11) increase. The increase 
of pi and p 2 is almost synchronous. Namely, as soon as t reaches the value (j — 1 -f aj)T for a particular 
integer j , a new term ipj gets included into the sum (2.11), i.e., the upper summation bound p> changes 
from its previous value j — 1 to the new value j. Similarly, as soon as t reaches the value (j — 1 + aj)T + Tint 
for a given jf, the term ipj drops from the sum (2.11), i.e., the lower summation bound changes from its 
previous value j to the new value j + 1. As has been mentioned in Section 2.3, in so doing the variation 
of the difference between p 2 and p\ never exceeds one. Moreover, the temporal interval that precedes the 
actual moment t and that is taken into consideration by formula (2.11) is again Tint- Consequently, we can 
still compute everything in the periodic framework, because the period Z (see (2.6)) is sufficiently large to 
accommodate the extent of retardation Ti nt , and as has also been mentioned it does not matter where on 
the period the computation of every given term starts. 

From the standpoint of implementation, when the upper bound p 2 increases by one at t = (j — 1 + aj)T 
nothing special needs to be done. If we simply continue marching equation (2.1) in the aforementioned 
periodic framework, then we will automatically start taking into account the new component of the RHS fj 
after t — (j — 1 4- &j)T. The situation with the lower bound p\ is somewhat different. Once it has increased 
by one (from j to j + 1) at t = (j — 1 + aj)T + T \ nt , the term tpj no longer needs to be included into the sum 
(2.11). However, in contradistinction to the case of Section 2.3 when all tpj were supposed to be computed 
independently of one another, here we cannot just stop the computation of a given (fj at t = (j - 1 +aj)T +T\ nt 
and subsequently say that <pj(x, t, T) = 0 for x € S(t) and for t > (j — \ + aj)T +T \ nt . Indeed, time-marching 
of equation (2.1) implies that all fragments of the solution < pj{x,t,T ) are calculated together as a sum and 
cannot be explicitly told apart. On the other hand, if we do nothing at t = (j — 1 + crj)T 4- T mt and continue 
with the time-marching, i.e., if we do not discontinue the computation of ^j(x, t 1 T) at t = (j — 1 +aj)T + T int 
and leave this term in the solution <p(x, t ), then right after this moment of time the first waves generated by 
fj at t = ( j — 1 + crj)T will start re-entering the domain 5(f) having traveled all the way across the auxiliary 
domain [E\ , E 2 ]. In other words, in the framework of the continuous time-marching with periodic boundary 
conditions, the term y?j(x,f, T ) cannot be left in the solution as it will “contaminate” the results on S(t). 

To avoid the aforementioned contamination, i.e., to prevent the re-entry of waves into 5(f), each term 
ifj(x,t,T) needs to be eliminated from the overall solution on the auxiliary domain [E \ , E 2 \ when the extent 
of its retardation (counted from inception) becomes exactly T, nt . For a given j the proper moment of time for 
elimination of ^(x,f,T) is t = (j - 1 + aj)T-h T int . Once we take out ^(x,f,T) at t = (j - 1 +oj)T + T int , 
this term may obviously be considered zero everywhere on [Ei^E^ for all subsequent moments as well. To 
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take out the term pj(x,t,T ) we need to interrupt the time-marching at t = ( j - 1 + aj)T + T\ ni , then go 
back to the inception moment of fj (x, t , T), which is t = (j - 1 + (rj)T } and independently integrate problem 
(2.10) for a particular j on [E\^E 2 ] from t 0 = (j — 1 + crj)T to £2 = 0 — 1 4- <rj)T 4- 7i nt . The result should 
then be subtracted from the time-marching solution at t = in the correct sense, i.e., both <p and ^ (rather 
their discrete counterparts) should be affected. Alternatively, we may notice that when integrating problem 
(2.10) from t 0 = ( j - 1 4- crj)T till t 2 — {j - 14- aj)T 4- Tj nt , the wave equation will, in fact, be homogeneous 
on a substantial portion of this time interval because fj(x,t,T) = 0 for t > t\ — (j+ 1 + aj)T. Consequently, 
instead of marching equation (2.10) over the entire time interval of length T\ nt , we may actually march it 
only from to = ( j — 1 + <rj)T till t\ = ( j 4* 1 +crj)T, then Fourier transform the discrete solution and advance 
it till t 2 = {j — 1 + crj)T 4- Tint by raising the corresponding amplification factors to the appropriate power. 
Numerically, this approach appears much cheaper, especially if it relies on highly efficient FFT subroutines. 

The new version of the lacunae-based algorithm has obviously been designed so that to exactly repro- 
duce the solution obtained with the original version of Section 2.3. The only difference is in the method 
of computation: Continuous time-marching in the periodic setup with cyclic subtractions of the retarded 
contributions versus separate computation of partial solutions driven by different components of the RHS. 
Consequently, the new version will posses the same properties as the original one. Foremost, it will provide 
for the temporally uniform grid convergence. Besides, it will obviously have linear computational complex- 
ity with respect to the grid dimension. (The cost of the FFT-based evolution in time distributed over the 
corresponding number of time steps is even less than linear if calculated per time step.) Finally, the algo- 
rithm will be universal in the sense that one will be able to build it as a modification of any consistent and 
stable finite-difference scheme. It will preserve the convergence rate of the original scheme while making the 
convergence uniform in time. 


2.5. Numerical Demonstrations. To actually demonstrate that the lacunae-based algorithm is an 
appropriate procedure that does deliver according to its theoretical design properties, we present some 
numerical results for the wave equation. For our simulations, we assume axial symmetry and employ the 
(r, z) cylindrical coordinates so that to account for the three-dimensional effects using two-dimensional 
geometry. Accordingly, equation (2.1) becomes: 


d 2 <p 
dt 2 


■<r - 


ld_ 

r dr 


dcp 

dr 


dz*J 


f{r,z,t) , t > 0 . 


(2.15) 


The solution <p of equation (2.15), as well as the RHS /, are functions of r, z, and t. The initial conditions 
for equation (2.15) remain homogeneous as before, see (2.2). 

We introduce the rectangular auxiliary domain [0 ,R] x [— Z/2,Z/2] of variables (r, 2 ), this domain is a 
specific realization of [Ex, E 2 ] shown in Figure 2.1. The boundary conditions are periodic with the period Z 
in the z direction, and zero Dirichlet at r = R: 


(2.16) 


<p(r,z± Z,t) = (p(r,z,t) , 

<p(i?, z y i) = 0 . 

The mathematical formulation of the problem obviously requires no boundary conditions at r = 0. However, 
for the purpose of subsequently building a discrete scheme (see below) we notice that the natural assumption 
of <p(r, z,t) being a bounded smooth function, along with the axial symmetry, immediately imply that 
= 0. Consequently, the Taylor expansion for p near r — 0 yields: 

1 <9 VI 


d<$ 

dr 


r=0 


vVv) = *(0,0 + 


2 dr 2 


• r 2 + 0(r 3 ) , 


r= 0 
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which means that 


dp 

dr 


d 2 p 
dr 2 


• r + 0(r 2 ) . 

r= 0 


Substituting the latter expression into (2.15) and considering the limit r — > 0, we obtain that on the 2 -axis, 
i.e., at r = 0, equation (2.15) reduces to: 


d 2 p 

dt 2 


o 9 2 >p 

dr 2 dz 2 


= , t = 0 , t > 0 . 


(2.17) 


To assess the quality of our numerical method we need to build a reference exact solution of problem 
(2.15), (2.2). This solution is obtained using the Lorentz’ transform: 

1 k/c z 


e = 


c = - 


k/c 


• t 


y/l - fc 2 /c 2 C ’ 
1 


(2.18) 


\/l — fc 2 /c 


* ct + 


\/nw 


Transformation (2.18) introduces the new coordinate system (r,£,0). The origin of this new coordinate 
system moves with the speed k along the 2 -axis of the original coordinate system. In other words, at every 
given t it is positioned at z — kt in the original frame of reference. In implementing transformation (2.18), 
we will always need to assume that k < c, as has also been suggested in Section 2.2. 

The key property of the Lorentz’ transform (2.18) is that it does not change the form of the wave 
equation (2.1) (and consequently, (2.15) and (2.17)), sec, e.g., [25]. As such, let us introduce an arbitrary 
function of time \ = x(^ x(0 = 0 for t < 0, so that it also be smooth on the entire R. Next, we define 
p 2 = r 2 + C 2 , and then 

X (9 - £ ) 

tf(r,C,*) = — (2.19a) 

P 


becomes a solution to the wave equation in the new coordinates (r, C#)- Solution (2.19a) is driven by a point 
5-type source, which is located at the origin {r = 0, £ = 0} and modulated in time by the function x{9). As 
X^O) = 0, this solution also satisfies the homogeneous initial conditions. Consequently, the function 


ip{r,z,t) 


p(r,z,t) 


(2.19b) 


obtained by substituting (2.18) into (2.19a) is a solution of equation (2.15) with the RHS /(r, z>t) = x(0 ’ 
5(r,z - kt). In other words, ^(r,2,£) of (2.19b) is a solution to the wave equation excited by a 5-source 
that moves straightforwardly and uniformly and is modulated in time by a given smooth function. Solution 
(2.19b) also satisfies homogeneous initial conditions (2.2). From the standpoint of physics, solution (2.19b) 
can be characterized as radiation of spherical waves by a moving point source. 

Solution (2.19b) is obviously singular. To use it for testing the numerical algorithm we need to remove 
the singularity. For that, let us first introduce the actual domain S(t). In all the experiments that follow 7 , 
the domain S(t ) is a sphere of diameter d with its center at the origin of the new 7 coordinate system: 
{r = 0, 2 = kt}. As such, this spherical domain moves uniformly along the 2 -axis, w 7 hich obviously helps us 
keep the axial symmetry intact. As has been mentioned, the speed of this motion is “subsonic,” k < c, which 
conforms to one of the key requirements for building the lacunae- based algorithm previously put forward in 
Section 2.2. 
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Let us also define r 2 = r 2 + (z — kt) 2 and introduce the function Q = Q(r ), f > 0, such that Q(0) = 0, 
Q(f) = 1 for f > Acrf/2, where k < 1, and also d = 0 for m = 1, 2, . . . till at least m — 4. 

Then, it is easy to see that the function ip{r,z,t) = ■ Q(r) is regular (continuous and bounded) 

everywhere. Moreover, it is easy to verify by direct differentiation that the same is true for the function 
f(r y z,t) 1= f □ </?(r, z,t), where □ denotes the wave operator, i.e., the left-hand side of equation (2.15). We 
will use /(r, t) defined this way as the source function for equation (2.15). Clearly, /(r, z, t) may, generally 
speaking, differ from zero only on the ball of a smaller diameter Kd concentric with S(t). Everywhere else, 
i.e., for f > Kd/ 2, /(r, z,t) = 0. 

Obviously, the solution of problem (2.15), (2.2) driven by this f{r,z,t) = D^(r, z,t) is the foregoing 

<p(r, Zj t) = z, t) • Q(r) . (2.20) 


This function satisfies the non-homogeneous wave equation with the RHS /(r, z, t) on a smaller ball of 
diameter Kd concentric with S(t). Everywhere else it is a solution to the homogeneous wave equation 
because it coincides with ip(r,z,t) of (2.19b). Consequently, <p(r, z,t) of (2.20) can be interpreted as the 
radiation of waves by a compactly supported moving source /(r, z,t). Numerically, we will be reproducing 
solution ip(r 9 z,t) given by (2.20) on the domain S(t) using finite-difference methods. 

We employ three different explicit central-difference schemes in our simulations. In all three cases we 
construct a uniform rectangular grid on the plane (r, z): r/ = lh r > l = 0,1,... , A T r , h r = R/N r , and 
z m = rnh z , m = 0, ± 1 , , ±A\, h z — Z/2N z . The discrete time levels are t n = nr, n = 0,1, — For the 
cell-centered second-order scheme, we keep the values of the unknown function p at the grid nodes in the z 
direction and at mid-points in the r direction: 


^r+1/2 ,m + ^/+ 1/2, 


n— 1 


1 1 


^/+3/2,m ^+l/2,m ^Z+l/2,m ^f-1/2,1 

n+i 7 n 7 

fir 


r / + l/2 ^ 

^«Vl/2,m+l “ ^H-l/2,m + ^f+l/2,m-l \ _ /n 

_ ~ ft 2 J " /i + 1/2 ‘ 


(2.21a) 


Equations (2.21a) hold for all / > 0. As in this case w^e do not have the unknown function defined on the 
axis of symmetry, and the closest values that correspond to l — 0 are half-grid-size away: p\ y 2jfn , then the 


scheme for l = 0 is obtained by simply assuming that ^±212^^ 


1=0 


= 0, which can be interpreted as 


a second-order approximation of the natural condition =0. This immediately yields for l = 0: 






1/2 ,m 


r— 0 


( 1 1 ^3/2, m Vi/2,m , ^?/2,m+l “ + ^1/2, m-I ^ _ /n 

WIC 1 K + V > - ' 


(2.21b) 


^ 7*1/2 ftr 

For the node-centered second-order scheme, is taken at the actual grid nodes, and for / > 0 we have: 


'P?£-W, m + 'Pl 


n — 1 


r/ h r 


r i+ 1/2 


V’f+l.m - 


•n-1/2- 




V?,m+1 - 2< P?,m + <P?,r> 

hi 


-) = /"« 


+ 


(2.22a) 
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To obtain the scheme on the axis of symmetry / = 0 in this case, we need to approximate equation (2.17). For 
the derivative in this equation we can first formally write I « ■ This expression 


obviously reduces to 
obtain: 


I r=0 


r= 0 


) VU.m — <^ 0,1 

- fcl! 


because of the symmetry: and consequently, we 


^m 1 2<p” m + 


2 A^ljTn ^Q,m + l ^^0,m — I A , n 


‘ 4 ft? ’ + ft? 


) = /o"r 


(2.22b) 


The last scheme is the node-centered fourth-order scheme. More precisely, it approximates spatial derivatives 
with the accuracy Oijtl + h *) and temporal derival ive with the accuracy 0{t 2 ). For l > 1 we have: 


¥>1 


n-f-1 




4 1 1 
3 T[ hr 
1 1 1 
3 r; 2 hr 


r /+ 1/2 


yi+l,rn ~ n,m 


h r 


r l- 1/2’ 




^/+2,m “ ~ VT- 2,m . . 

n+i l xt n-i 1 + 


2h r " 1 2h 

-y’tm+Z + 16 ^m+l - 3(V” m + 16^,! - y” ;Wi _ 2 
12/if 

For l — 1 we have r/_i = ro = 0 and consequently: 

¥>"£ - M" + j 


(2.23a) 


)=/"» 


/ 4 1 1 

{ Pl > m ~ ‘Po.m 

11 1 L <?3, m -^"ml 

\3 r\ h r 

[ r3/2 hr Fl/2 hr J 

3 ri 2ft r [ 2 h r 


-^"m+2 + 16^” m+ l - 30<p? m + 16< m _ x - V3" m _ 2 

12/i? 


) = /m 


(2.23b) 


Finally, for l = 0 we again have to approximate equation (2,17). Using symmetry like in the previous case, 
we arrive at: 


«-2¥>o" m +^,m 




T Z 

+ 32^",; 


30<^o . 


12/i2 

-yo,m+2 + IggjU+l - 3(Vg m + 16yp m _] - ¥?g [fn _ 

12/i? 


-) = / 0 "r 


(2.23c) 


For all three schemes, (2.21), (2.22), and (2.23), setting the discrete boundary conditions (2.16) on the outer 
boundary of the auxiliary domain [0,R] x [— Z/2,Z/2] is straightforward. An extra boundary condition is 
needed for the fourth-order approximation. As it basically does not matter what boundary conditions we 
use on the outer boundary of the auxiliary domain (see Section 2.2), we simply set m = 0 in addition 

to <fx r m — 0. Regarding the time step r, all three schemes are explicit and as such, there is a Courant-type 
stability constraint. 

As has been mentioned, we present the results of numerical computations that follow' in order to corrob- 
orate the theoretical design properties of the lacunae-based algorithm, i.e., to show' the temporally uniform 
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grid convergence on long time intervals. For that purpose, we conduct a grid refinement study, i.e., approxi- 
mate the exact solution (2.20) on a sequence of successively more fine grids. In so doing, the time step r for 
the two second-order schemes (2.21) and (2.22) is always reduced with the same rate as the corresponding 
spatial sizes h r and h z \ and the time step r for the fourth-order scheme (2.23) is reduced twice as fast 
(i.e., by a factor of four every time h r and h z are reduced by a factor of two) so that to demonstrate the 
fourth-order overall convergence in the end. The computations in each case were run till the dimensionless 
time t = 200 * d/c, i.e., for 200 times the time interval required for a w r ave to cross the domain. This certainly 
qualifies as “long term” from the standpoint of any conceivable application. 


Lacunae-based (2,2) cell-centered scheme Lacunae-based (2,2) node-centered scheme 



Dimensionless time Dimensionless time 


(a) The second-order scheme (2.21) (b) The second-order scheme (2.22) 

Fig. 2.2. Grid convergence study for the long-term lacunae-based integration of the wave equation. 


Lacunae-based (2,4) node-centered scheme 



Dimensionless time 

Fig. 2.3. Same as Fig. 2.2 for the fourth-order scheme (2.23). 


The actual parameters that we have used for 
our simulations are the following: R — Z = 7r, 
d = 1.8, c = 1, k = 0.2, k = 0.8. The spatial grid 
is composed of square cells: N r = N z and conse- 
quently, h r = h z = h. The actual grid dimensions 
N r x 2N~ are: 64 x 128, 128 x 256, and 256 x 512. 
The temporal partition size 2T, see (2.8), is found 
from (2.6) assuming that T\ n t = ; the 

overlap parameter a — 1/2. The functions 0(t) 
and Q(r) on the intervals of their variation from 
0 to 1 are built as polynomials of degree 9 (with 
only odd powers included), which guarantees four 
continuous derivatives in transition to the constant 
(either 0 or 1). The function \{t) is defined as fol- 
lows: x(f) = (l + | sint) P (1 — 2 ~) , wdiere P(t) 
is, again, a polynomial of degree 9 that decays 
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smoothly from 1 to 0 on the interval [0, 1] (four continuous derivatives). Finally, to subtract every <pj 
from the overall solution at the proper moment t = (j — 1 -f crj)T + T int , we first march equation (2.10) 
from t — (j — 1 + a j)T till t = (j + 1 + oj)T and then use Fourier expansion in z and expansion with 
respect to the corresponding discrete eigenfunctions (calculated numerically) in r to advance it further till 
t = {j — 1 + &j)T + Fjnt* 

In Figure 2.2 we show error profiles (more precisely, natural logarithm of the relative error on the domain 
S(t) in the maximum norm as it depends on the dimensionless time) on all three grids for both second-order 
schemes (2.21) and (2.22). In Figure 2.3, similar curves are shown for the fourth-order scheme (2.23). From 
these figures we conclude that indeed no error is accumulated in the course of computations because all error 
profiles are flat throughout the entire 200 * djc time interval. Thus, the solution does not deteriorate as 
time elapses. Figure 2.2 also shows that every time the grid is refined by a factor of two the error drops by 
approximately a factor of four, which indicates the second-order convergence. Similarly, Figure 2.3 shows that 
every time the grid is refined by a factor of two the error drops by approximately a factor of sixteen, which 
is an indication of the fourth-order convergence. Consequently, we can conclude that numerical experiments 
fully corroborate the theoretical design properties of the algorithm. 

3. Lacunae-Based ABCs for 
the Wave Equation. The lacunae- 
based algorithm of Section 2 provides 
a venue for constructing the ABCs for 
a class of problems that reduce to the 
homogeneous wave equation in the far 
field. We schematically depict the ge- 
ometric setup for one such problem 
in Figure 3.1, assuming for simplic- 
ity that there is no source motion, 
k = 0, and the computational domain 
is stationary. We emphasize though 
that this is not a limitation, and that 
the actual ABCs will be constructed 
and tested for the general case of a 
moving computational domain, w T hile 
the law 7 of motion can be arbitrary, 
see Section 2.1. The problem to be 

solved on the bounded interior do- 
Fig. 3.1. Schematic geometric setup for the ABCs. maillj i. e ., in the near field, see Fig- 

ure 3.1, may involve some complex 
phenomena whose nature, however, is not essential for the current discussion. 4 We only require that the 
overall combined formulation of the problem be uniquely solvable and well posed under the assumption of 
radiation of waves in the far field (from S(t) tow r ard infinity), where the problem is assumed to be governed 
by the homogeneous wave equation. The role of the ABCs (as mentioned in Section 1) is to provide a closure 
for the truncated problem solved on the actual computational domain S(t). This closure has to ensure that 

4 The interior domain is, of course, the same as S(t) of Section 2; for the stationary case we obviously have S(<) = S(0). 
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the corresponding finite-domain solution recovered with the help of the ABCs be close to (ideally, exactly 
the same as) the solution of the original non-truncated problem restricted to the bounded domain, see [2]. 

In the wave propagation framework adopted in this paper, one can say that the ABCs have to replace 
the entire far field, i.e., everything beyond the bounded interior domain S(t), so that the resulting artificial 
boundary be completely transparent for all the outgoing (i.e., radiated) waves. We also note that the 
incoming waves, provided that they are meaningful for a particular setup, can, in fact, be taken into account 
through the boundary conditions as well, but we do not discuss this issue here for the reason of simplicity. 

3.1. Preliminary Considerations in the Continuous Framework. Let tp c = <Pc{x, t) be a solution 
to the aforementioned combined problem. In the far field, i.e., outside 5(t), the function <p c {x, t) satisfies the 
homogeneous wave equation. We also assume for simplicity that the solution <p c (x,t) “smoothly originates 
from zero” at t = 0, (i.e., turns into zero along with its first derivative) in much the same w T ay as the 
solution </?(#,£ ) of (2.1), (2.2) does. This assumption, in fact, will present no limitation when constructing 
the ABCs. The argument is the same as the one that allows us to relax the assumption of homogeneity of 
initial conditions when building the original lacunae- based algorithm, see [23]. 

Let us now introduce a special multiplier function that is again schematically shown in Figure 3.1. This 
function /r = fi(x^t) is defined for all those x and £, for which the solution ip c (x,t) makes sense. We first 
require that > 0,Vx £ S(t) : v(x,t) = 1, or in other words, that the multiplier be identically equal to 
one everywhere outside the computational domain S(t) for all times. We also require that the multiplier 
be identically equal to zero, /i(x,f) = 0, on most of the domain S(t) (again, for every t) except next to 
its boundary from the interior side. An example of the narrow near-boundary transition region, where 
the multiplier fi{x t t) changes its value from zero to one, is shaded on Figure 3.1. What is important, we 
require that the multiplier ^ i(x,t ) be a sufficiently smooth function with respect to both x and t , which 
essentially means that the transition within the shaded region on Figure 3.1 has to be smooth. Regarding 
the time dependency of n(x, £), once the domain S(t) moves according to a prescribed law, the construction 
of the multiplier has to trace that motion. If the computational domain is stationary, S(t) = 5(0), then the 
multiplier still may, but does not have to, depend on time. 

Next, we apply the wave operator □ = Jp- — c 2 A of (2.1) to the function n(x,t) - <£> c (:r,£), which is 
defined everywhere, i.e., both inside and outside S(t ). 5 We will obviously have: 


= 


d 2 fiip c 

at 2 


— c 2 A(/i<£ c ) = g{x,t) 


r 

= 0 

= 0 


Vi, Vx # S{t) 

in the transition region 

“well inside” S(t) 


(3.1) 


The function g(x,t) of (3.1) may generally speaking differ from zero only in the foregoing near-boundary 
transition region; it is zero outside S(t) because the function f.up c coincides there with the solution ip c of 
the homogeneous wave equation; it is also zero inside S(t) because ji~ 0 there. On Figure 3.1 the non-zero 
portion of g{x i t) is identified as the right-hand side RHS. 

We can now consider the problem (2.1), (2.2) with the function g{x,t) of (3.1) substituted instead of the 
generic RHS f(x, t). The key fact that we will need for constructing the ABCs, and that follows immediately 
from the unique solvability of the Cauchy problem for the wave equation, is that the solution to this problem 
will coincide with n(x,t) • <p c {x,t ) everywhere. What will be of particular importance to us is that as such, 


5 Note, the solution <p c (x,t) may not be defined on all of S(t) if, e.g., there is a scatterer inside. As, however, = 0 

there, we can consider t) to be defined everywhere. 
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this solution will coincide with ip c (x,t) outside S(t) for all times, because p(x,t) = 1 there. In other words, 
we have replaced all of the interior problem on S{t ) (no matter how complex it may be) by the special 
near-boundary source function g(x,t) so that the corresponding far-field portion of the solution, i.e., the 
wave-radiation solution outside S(t ), remains totally unaffected. Later on, see Sections 3.2 and 3.3, this 
reduction interpreted in the discrete framework will be used for setting the ABCs. The idea is to use the 
exterior solution obtained in an alternative way through integrating the near-boundary sources as a closure 
for the interior problem solved on the finite computational domain. 


3.2. The Concept of Discrete ABCs. To construct the ABCs for a finite-difference scheme that 
approximates the problem described in the beginning of Section 3, we will employ the considerations similar 
to those of Section 3.1, but on the discrete level. As a helpful illustration, we will first consider here a 
one-dimensional model example, 6 and then, in the next Section 3.3, show how to build the ABCs for the 
actual multidimensional wave propagation problems. 
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Assume that we are solving a one- 
dimensional combined problem on the entire 
R The computational domain S(t) = S (i.e., 
the “near field” ) is fixed, it is the half-line x < 
0 (more precisely, it is {(#, £)| x < 0, t > 0}); 
its complement M \S represents the “far field,” 
which is to be truncated and replaced with the 
ABCs. As such, the ABCs are to be set at the 
interface x = 0. In accordance with the pre- 
vious discussion, we also assume that the far 
field is governed by the one-dimensional ho- 
mogeneous wave equation 


d 2 tp 2 d 2 <p 

w~ c ox* 


= o, 


(3.2a) 


which is approximated by the standard 
second-order central- difference scheme 


Fig. 3.2. Illustration to the one-dimensional example. 






,2 fj+l ~ 2^7 + Vj-l 

h 2 ’ 


(3.2b) 


constructed on the rectangular grid of variables x and t with sizes h and r = hjc respectively, using the 
five-node stencil shown in Figure 3.2. (Note, all the schemes used for simulations in Section 2.5, see (2.21), 
(2.22), and (2.23), are of the same central-difference explicit three-level type. This, however, is by no means 
a limitation — the ABCs can be constructed for any type of discretization.) 

To create the discrete near-boundary sources similar to those of Section 3.1, and eventually set the 
discrete ABCs, we will need to be able to apply inside S the same finite-difference wave operator of (3.2b) 
as the one we are using on R\S. As such, we formally extend the exterior discretization, i.e., the rectangular 
grid with h x r cells, into the interior domain 5, as shown in Figure 3.2. We re-emphasize, however, that 


6 Generally speaking, one-dimensional problems do not have lacunae (except in special cases); as such, this example will 
only demonstrate the formal construction of the ABCs on the grid. 
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this is done only for the “artificial” purpose of building the ABCs. The actual governing equation in the 
near field, i.e., on 5, as well as its discrete counterpart, may be more complex than the wave equations (3.2) 
above. In fact, neither the scheme stencil nor the grid used for computations in the near field have to be 
the same as those in the far field (although, they, of course, may). We only require that the exterior scheme 
(3.2b) with the stencil shown on Figure 3.2 be applicable till x = 0. And all we need to know from the 
standpoint of setting the ABCs, is that the overall combined problem (near field and far field) is uniquely 
solvable and well-posed. 

Let us now consider all nodes of the aforementioned rectangular grid that belong to 5 (on all time levels), 
i.e., those for which xj = jh < 0. We denote this set of grid nodes by A r+ ; the complementary set that 
consists of all nodes from E \S, i.e., those, for which xj > 0, is denoted by A r ~ , see Figure 3.2. If we formally 
apply the five-node stencil of the scheme (3.2b) to each node from A/’ + , then this stencil is obviously going 
to sweep one more vertical row of nodes, which already belongs to (i.e., to R\S), and which is denoted 
by 7“ on Figure 3.2. Reciprocally, if we apply the stencil to every node in A/"”, it will also sweep the nodes 
7 + that are already in A r+ . The two-layer grid structure 7 = 7+ U 7“ will be called the grid boundary; it 
represents on the discrete level the continuous interface between S and K\S, which is the vertical line x = 0. 

Next, w T e assume that we integrate the interior problem one time step after another, and that we already 
know the discrete solution on the domain 5, as well as the values of on the grid boundary 7, up to 
a certain time level n (in particular, n may be equal to zero, which corresponds to the initial conditions). 7 
These data obviously allow us to advance the next time step n + 1 on 7 + and everywhere inside 5; in so 
doing, the outermost interior location 7+ on the level n + 1 is computed by scheme (3.2b) using the stencil 
shown in Figure 3.2. These data, however, already do not allow us to calculate the discrete solution ipW 
at 7” on the level n - hi. And if this value p™- 1 is not knowm, then we cannot advance further to level 
n -{- 2. Therefore, we conclude that the function of the ABCs in the discrete framework will be to provide 
the missing boundary values of the solution at 7“ on all time levels, one after another, starting from n = 1. 
This indeed constitutes the closure of the discrete system solved on 5. 

To provide the foregoing missing boundary value for a given n, we recall that even so we do not 
know the discrete solution on level n + 1 beyond 7+ (i.e, we do not know ^t 1 ), we do know that the 
solution can be complemented on A r ~ to a solution of equation (3.2b) on all time levels till n - j- 1. For 
our purposes, w f e will only need the existence of this complement rather than its actual representation. Let 
us now introduce a multiplier function p similar to that we have used in Section 3.1. The near-boundary 
interior transition region for this multiplier is schematically shown by the shaded area on Figure 3.2. We 
apply this multiplier to the combined discrete solution = p^\ UJyj r_ on all time levels including n + 1 
(obtaining may require projecting onto the grid A r+ ). In so doing, nothing changes on J\f~ U7 + , 
because /i = 1 for x > 0. All the changes due to multiplication of by p will obviously be introduced 
on A r+ \7 + only. Those amount to a smooth passage within the transition region (see Figure 3.2) from the 
actual unaltered values of the solution on 7 + to zero “well inside” the computational domain. 

Next, similarly to Section 3.1, we apply the discrete wave operator IlK^ of (3.2b) to the modified solution 
fi(fc h \ As is defined up to the level n + 1, the result will be defined up to the level n. Analogously 


7 We use the subscript “5” 
on a different grid. 


in rather than ” to emphasize that the actual interior discrete solution may be computed 
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to (3.1), we obtain for all levels till n: 

r 

= 0 on the grid A r ~ 

^ 0 on 7 + & in the transition region (3.3) 

= 0 “well inside” A r+ 

Notice, we can claim that the result in (3.3) is zero on M~ only because nothing has been modified by ft on 
7 + and beyond. As such, we are simply using the fact that is a solution to the homogeneous equation 
outside the computational domain, and we do not need to know the explicit form of this solution. 

Finally, we solve the non-homogeneous counterpart to equation (3.2b) driven by the RHS of (3.3) 
everywhere on A r+ UA r “; this will be henceforth referred to as solving the auxiliary problem. According 
to our construction, solving the auxiliary problem will allow us, in particular, to recover the value of 
which was not known previously, and which can now be supplied to the interior scheme as the missing 
boundary value. This means that we will have provided the ABCs for the interior problem, because in so 
doing we complete the time level n + 1 and facilitate advancing the next level n + 2. 

The are several important comments to be made regarding the foregoing ABCs’ algorithm. At a first 
glance, the new formulation simply does not change much from the standpoint of solving the original infinite- 
domain problem. Indeed, all we have done is replaced the interior problem by the artificial near-boundary 
sources so that the exterior solution remain unaffected. Then, we suggested to use this exterior solution 
to close the interior discretization. However, obtaining this solution, i.e., solving the auxiliary problem, 
basically brings along the exact same set of complications that we have been trying to avoid by introducing 
the ABCs. Indeed, the domain of the auxiliary problem A r+ UAf - is still unbounded and as such, special 
treatment will be required for its numerical solution. 

There is, however, a fundamental difference. The new auxiliary problem is linear throughout the entire 
space, and it is driven by known sources that are compactly supported inside the computational domain S. 
Consequently, the lacunae-based algorithm of Section 2 appears to be a most natural tool to solve it. 8 Em- 
ploying the lacunae-based algorithm immediately implies that the domain of the auxiliary problem becomes 
bounded. Moreover, the “sufficiently retarded” sources do not contribute to its solution (see Section 2), i.e., 
only a limited extent of temporal pre-history of the solution will be needed to sustain the continuous time 
marching no matter how far in time we would like to advance the solution. In other w r ords, the missing 
boundary value for the interior discretization tp^t 1 can be obtained using only finite computer resources in 
terms of both memory and number of arithmetic operations. Furthermore, these resources (say, per time 
level) will not increase no matter for how long we may need to run the computation, i.e., how large n may 
become. In this sense, the proposed ABCs become “true ABCs,” i.e., the procedure that guarantees the 
appropriate closure of the truncated problem with only finite non-growing amount of computer resources 
required. In addition to that, we are guaranteed that the ABCs as a part of the overall algorithm will not 
contribute toward the buildup of numerical error during long runs. 

The proposed ABCs can obviously be implemented via alternating interior /exterior steps. Namely, we 
advance one time step in the interior (including 7 + ) assuming that all the data that we need from the 
previous time levels are available. The resulting newly calculated time level will be the only one to which the 
multiplier has not been applied yet. We multiply it by fi 7 and then apply the direct operator thus obtaining 

8 We reiterate that this algorithm cannot be applied in the case of one space dimension, but our ultimate goal is three- 
dimensional problems, see Section 3.3, and the considerations of the current section are for illustrative purposes only. 
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the right-hand side g^ h \ see (3.3), on one more time level as well. Finally, we perform one step of the 
lacunae-based integration of the auxiliary problem driven by g ^ and obtain the missing boundary value for 
the interior problem. Then, the procedure cyclically repeats itself. Summarizing, we can say that having 
advanced the interior solution, we then generate a new contribution to the RHS of the auxiliary problem 
and subsequently advance its solution, which, in turn, allows us to calculate the next interior step. 

An important observation, which is easy to make, is that the missing boundary value that the 

ABCs provide does not, of course, depend on the actual shape of the multiplier g in the transition region. 
Indeed, g is defined so that it does not alter the solution on the grid boundary 7 . Consequently, when we 
first apply the direct operator to gp>[ h \ see (3.3), and then integrate the non-homogeneous wave equation 
driven by g^ h \ the solution on 7 will remain unchanged no matter what changes have been introduced by g 
in the interior. As such, the value y?"t* will only depend on the values of the solution on 7 on all previous 
time levels, as well as on y?”* 1 . Moreover, since all the operations that we perform when constructing the 
ABCs are linear, we can symbolically write the resulting boundary condition as a linear form: 


<P 


71+1 

7 “ 


= l(v> 


71+1 
7 + 




71—1 
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(3.4) 


Technically, the dependency of on the previous time levels, see (3.4), involves all of the latter, from 
ip™ all the way back till n = 0. However, the use of the lacunae in three space dimensions will allow us 
to truncate (3.4) and leave only several levels that immediately precede n + 1; the number of the levels 
involved will be fixed and will not increase with the increase of n. As such, temporal non-locality of the 
ABCs will be limited, and this will not be a consequence of any approximation, but rather an implication of 
the fundamental properties of the problem. We also note that the representation of the ABCs in the form 
(3.4) serves primarily the reason of convenience and compactness in notations. In fact, the coefficients of 
the linear form l never need to be known explicitly except, possibly, when multiple interior problems are 
solved with the same exterior model (means same grid, same geometry, same scheme). In this case it may 
be beneficial to calculate the form l once and ahead of time, compared to the straightforward calculation of 
many times according to the procedure outlined above. 

It is also important to mention that boundary condition (3.4) can be obtained in the framework of a 
general unsteady ABCs’ methodology proposed by Ryaben’kii in [24] (see also older work [26]) for a variety 
of problems, including multidimensional cases, domains of varying shape, and different types of schemes — 
explicit as well as implicit. Work [24,26] describes the theoretical construction of the ABCs per se, and does 
not address any issues related to the actual computations (for example, using lacunae-based integration, a s 
proposed in the current paper). The methodology of [24,26] relies on the concepts of generalized potentials 
and boundary projection operators of Calderon’s type obtained and implemented in the discrete framework 
by means of the difference potentials method, see [27-29]. In this perspective, the ABCs of [24,26], and 
boundary condition (3.4) in particular, can be interpreted as discrete counterparts to Calderon’s boundary 
equations with projections in the unsteady case. In the following Section 3.3 we will describe a direct 
approach to obtaining multidimensional ABCs on moving boundaries, with no explicit use of the apparatus 
of Calderon’s projections, and will also show how r to apply the lacunae-based algorithm to perform the 
computations needed for these boundary conditions. 

To conclude this section, w r e emphasize that even so y^t 1 obtained according to (3.4) formally does not 
depend on the shape of the multiplier g inside 5, w r e still need to have this multiplier smooth. In other 
words, we could not have used, e.g., a step function, instead of g. The reason is that non-smoothness will 
ruin lacunae in the discrete solution (see Section 2 and [23] for more detail) and consequently, we will no 
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longer be able to use the lacunae-based integration for solving the auxiliary problem and as such, setting 
the ABCs. 


3.3. General Construction of Discrete ABCs. Similarly to the setup of Section 2, we now consider 
a domain S(t) € M 3 that has finite diameter d for all times, and other than that may travel in space according 
to a prescribed law with the only limitation that its maximum speed be subsonic: k < c. S(t) will be the 
computational domain, or near field. In the far field, i.e., outside S(t ), we assume that our model is governed 
by the homogeneous w r ave equation: 


d 2 <p 
8t 2 


V dxj dx% dxj ) 


= 0 , 


t > 0 . 


(3.5) 


As we have discussed, inside S(t) the solution <p — p{x,i) may be governed by a more complex equa- 
tion/system, but all we need to assume is that the overall problem be uniquely solvable and well posed under 
the condition of waves’ radiation toward infinity. For simplicity, we also assume homogeneity of the initial 
data everywhere, which is, however, not a limitation (see [23]). 

Let us now introduce the discretization grid for equation (3.5). In principle, we need this grid only in the 
far field, i.e., outside 5(t), because the interior problem may be discretized in a different way, as indicated 
before. As, however, we have also seen, to obtain the ABCs we need to set up the auxiliary problem for the 
non-homogeneous counterpart of equation (3.5) driven by the special near-boundary sources. The auxiliary 
problem is to be formulated and solved on the entire space. As such, we introduce the grid for the linear wave 
equation on the entire R 3 x [0, oo) as well. We denote by M the collection of all grid nodes in R 3 x [0,oc), 
on which we evaluate the solution ip. Since (3.5) is an evolution equation, it is convenient to consider A f as 

a composition of spatially aligned grid hyper-planes: Jsf — A/o U Aq U . . . U M n Each J\f n is a spatial grid 

on K 3 , and we emphasize that they may, but do not have to, be the same on different levels n. 

Let the individual nodes of the grid A/* be denoted by n. Equation (3.5) is approximated by a finite- 
difference scheme, which we assume, of course, to be consistent and stable: 


yi ClmnVn = 0. (3.6) 

neA/'m 

In (3.6), A r m denotes the stencil attributed to the node m, and a mn are the corresponding coefficients. When 
we say that the stencil is attributed to a particular node, we mean that the residuals of the discrete equation 
are evaluated at this particular grid location. In regard to this, we note that the residuals of the discretized 
equation (3.5) may, but do not have to, be evaluated on the same grid Jsf. To preserve the generality of 
the discussion, we assume that there is another, different, grid Ad in M 3 x [0, oo), on which we keep the 
residuals, as well as the right-hand sides, if any, of the discrete wave equation. The subscript m in equation 
(3.6) basically refers to this grid: m £ M. In the one-dimensional example of Section 3.2, both grids A f and 
M w r ere simply the same, and we did not have to distinguish between the two. To give an example of the 
opposite type, we mention the Yee scheme, see [30], which is one of the primary tools for discretizing the 
Maxwell’s equations, 9 and which involves staggering in both space and time. 

Next, w^e introduce two subsets of nodes of the grid Jsf. Let the level Af n correspond to the actual moment 
of time t n . For every n, we define Jsf+ as the set of all those and only nodes on this level that belong to the 
domain S(t n ), and A r “ as the complement of to the entire A r n : A"” = A r „ \ A^ . In other words, Jsf~ 

9 The simplest version of the Maxwell’s equations describes the propagation of electromagnetic waves in vacuum, which is 
a wave model similar in many respects to that given by (3.5). 
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contains all those and only those nodes of Af n that belong to M 3 \S(i n ). Subsequently, we define the set 
as the composition of all Af+ for all levels, and the set A r ~ as the composition of all for all levels: 

u k, j*-= u a;-. 


Clearly, = A r \A' + . 

In our definition of the scheme, see (3.6), we have identified the stencil M m with the grid location m, 
at which the residual is evaluated. Prom here on, we will assume for simplicity that the scheme (3.6) is 
explicit. In this case, there is only one non-zero coefficient a mn on the upper time level of the stencil j\T m . 
We will denote the corresponding grid node by n, and when it may not cause confusion, we will refer to the 
same stencil as either M m or Jsfh- If will also be convenient to introduce the four-dimensional (space-time) 
vector b = fi - m. The meaning of this vector is that is defines the relative position of the node n, at which 
the upper-level coefficient is non-zero: a m „ / 0, with respect to the “center” of the stencil m. This vector 
is obviously constant, it depends only on the local structure of the stencil, and does not depend on where 
exactly on the grid this stencil is applied at every given moment. In the one-dimensional example considered 
previously, we would have b = (r, 0). 

Let us now apply the stencil Mn to 
every node n € in so doing, the 
stencil obviously sweeps the entire grid 
, as w r ell as a portion of the grid Af~ 
next to the interface, we will denote this 
portion by 7“: 
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On Figure 3.3, we present a one- 
dimensional illustration 10 similar to 
that in Figure 3.2 for the case of a uni- 
form motion of the computational do- 
main, when the space-time trajectory of 
the boundary is a straight line (the set 
7“ is denoted by small circles). From 
the standpoint of implementation, the 
values of the solution at the nodes 7' 
are exactly those that need to be provided by the ABCs from the exterior side so that to be able to calculate 
the solution at every interior node n 6 A f + using the scheme (3.6). Reciprocally, the stencil A/f» applied to 
every node n E J\f~ sweeps additional nodes 7 + C A" + , see “bullets” on Figure 3.3: 


Fig. 3.3. One- dimensional illustration to the case of a moving domain . 


7 + = 


U #n 

nC.V' - 


fK 


The set 7 + complements 7 to the complete grid structure known as the grid boundary 7 (see [28]): 

7 = 7 + 1J 7~- 


10 We note again that everywhere in this section the one-dimensional examples are for illustration purposes only, the actual 
algorithm is three-dimensional. 
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The grid boundary is a multi-layer fringe of nodes (two-layer fringe in the particular case of a second-order 
scheme with the stencil depicted in Figure 3.3) that is located near the continuous boundary and straddles 
it in some sense (cf. to Section 3.2). 

To proceed with the construction of the ABCs, we will need to assume hereafter that the region of 
linearity, i.e., the area where the solution of the overall combined problem is governed by the linear homoge- 
neous wave equation, extends “a little bit” to the interior of the computational domain S(t) as well. More 
precisely, this region will be assumed to extend inward at least as far as the entire grid boundary 7. This 
obviously presents no limitation from any standpoint. The multiplier fi = /i(a:,f) in this case is required to 
be identically equal to one, /i(x, t) = 1, not only outside 5(f), but also inside — again, to the extent of 7+. 
As such, the transition region for the multiplier, which is schematically shown by the darker gray shading 
on Figure 3.3, is shifted away from the boundary of 5(f). 

To actually build the ABCs, we will perform the same procedure as outlined in Section 3.2. What we 
actually need in the discrete framework is to obtain the missing exterior boundary values of the solution 
on every time level n 4- 1, or in other words, to complete this time level, so that to be able to advance 
the next time step. To do that, we take the solution already computed inside 5(f) up to the level n + 1, 
multiply it by f.i and then apply the discrete operator everywhere. In doing so we assume, as mentioned 
before, that starting from 7 outward, the solution satisfies the discrete homogeneous wave equation. In the 
general case that we are looking at now, an application of the operator brings us from the grid A r to 
the grid M . The construction of the grid boundary 7 and multiplier \i guarantees that on all time levels up 
to n the near-boundary artificial sources g ^ will satisfy (cf. to (3.3)): 



= 9 (h) 


= 0 
^0 
= 0 


for such m € M that m + b € 

for such near-boundary m G M that m + b € A r+ 

one the grid M "well inside” 5(f) 


(3.7) 


Formula (3.7) suggests that in addition to the actual interface between the domains, i.e., the boundary of 
5(f), it will also be convenient to consider another space-time trajectory obtained from this original interface 
by the constant displacement —6, it is shown by the dash-dotted line in Figure 3.3. The righ t-hand side g ^ 
will be zero on M everywhere outside this new displaced boundary, and will differ from zero right next to it 
on the interior side. On Figure 3.3, we schematically show by the lighter gray shading the region where we 
still have /i = 1 but the RHS g ^ may already differ from zero. Note, in the example of Section 3.2 be did 
not need to consider the displaced interface because the domain 5(f) was stationary and the displacement 
b was parallel to the time axis: 6 = (r, 0). Another important case when the two boundaries would 
coincide is n = m & b — 0. However, we cannot generally assume that for time-dependent problems. 
On the other hand, we mention that the grid boundaries originally introduced in [27, 28] and previous 
publications by Ryaben’kii for the solution of steady-state problems using difference potentials method, have 
been constructed so that just the center m of the stencil (where the residuals are evaluated) would sweep 
a given grid subdomain and as such, generate the aforementioned fringe of nodes 7 next to the continuous 
boundary. 

Let us now make a few remarks of explanatory nature regarding the structure of the grid boundary 7. 
It obviously depends only on the type of the stencil Af m and geometry of the actual continuous boundary 
that it straddles. From the definition of 7 is is easy to see that once we have a solution to the homogeneous 
equation on 7 and everywhere in the exterior, and operate by on this solution, then we can guarantee 
without actual calculation that the result will be zero for all m G M: m + b e J\f ~ . However, we cannot 
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“touch” even one single node from 7 so that not to lose this property. If, for example, we allows for an 
alteration (via ft / 1) of a node from 7+ (see Figure 3.3), it will necessarily affect the exterior RHS. The 
latter may, generally speaking, become non-zero at some node(s) m E M: m + b E A r ~, and we will no 
longer be able to actually calculate it because we do not know anything about the exterior solution beyond 7 
except that it satisfies the homogeneous equation. We see, therefore, that it would ruin the entire derivation. 
On the other hand, the construction of the grid boundary 7 is consistent in the sense that to calculate the 
actual non-zero near boundary sources for those m E M, for which m + b E A' 4 *, it is sufficient to know 7 
the solution only on 7 and further inw T ard, nothing outside 7 needs to be knowm. As for the values that 
are still needed, those are provided by the ABCs’ algorithm on every time level and as such, are available 
on all subsequent levels for calculating the source terms of the auxiliary problem. 

Having outlined the construction of the grid boundary 7 and near-boundary sources g ^ in the general 
case, we build the actual ABCs’ algorithm in much the same w T ay as described previously. We perform the 
alternating interior /exterior steps: First advance one step in the interior, then apply /x and calculate one 
more level of the sources g^ h \ and finally make one step of the lacunae-based integration of the auxiliary 
problem driven by these sources (see Section 2), thereby providing the missing data for advancing the next 
interior step. Then, the procedure cyclically repeats itself. To solve the auxiliary problem in this general case, 
w T e will obviously need a full-fledged version of the lacunae-based algorithm (see Section 2) that accounts for 
the motion of the sources and employs periodic boundary conditions in space and continuous time-marching 
with cyclic subtraction of the retarded contributions. As has been shown, implementation of the lacunae- 
based integration technique guarantees that the domain of the auxiliary problem wall be bounded, and the 
computer resources needed for the ABCs will be finite and will not grow with time. Other properties of 
the ABCs outlined in Section 3.2, namely, independency on the shape of the multiplier /x, and possibility to 
express the boundary values on the current time level as a linear function of the values on the previous levels, 
see (3.4), hold in the general framework of this section as well. Because of the lacunae, the aforementioned 
linear form will depend only on the finite non-increasing number of the preceding time levels n, essentially 
those included in the summation (2.11) once this formula is discretized on the grid A r . This means that the 
temporal nonlocality of the ABCs wall be limited. As for the multiplier /x, it ha s to be chosen sufficiently 
smooth so that to maintain good quality of the lacunae in the discrete solution, see Section 4. 

4. Numerical Experiments with the ABCs. 

4.1. The Wave Equation with Known Exact Solution. The first case that we analyze in the 
framework of the ABCs is actually the exact same problem that we solved in Section 2.5. It is the wave 
equation driven by a compactly supported oscillatory source in straightforward uniform motion. The exact 
solution of the problem w r as obtained using the Lorentz transform. The key difference between the current 
approach and that of Section 2.5 is that previously we have applied the lacunae-based algorithm directly to 
the original problem. Here, we rather decompose the problem into the near field S(t) and far field M 3 \5(t), 
even so both are governed by the same wave equation. The integration in the near field is then performed 
by the conventional time marching. The exterior closure needed to sustain this time marching is provided 
by the discrete ABCs on the boundary of S(£). The ABCs are constructed on the basis of the procedure 
outlined above — through the lacunae-based integration of the artificial near-boundary sources. 

Similarly to Section 2.5, we have implemented three different schemes (2.21), (2.22), and (2.23), and 
every time integrated the problem till the dimensionless time has reached t = 200 • d/c . The multiplier 
fi = fi(r,z,t) was constructed so that to have four continuous derivatives w T it,h respect to all its arguments. 
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Smooth transition from 0 to 1 was obtained with the help of algebraic polynomials of degree 9 similar to 
those used in Section 2.5. The extent of the transition region has varied slightly between different cases with 
no noticeable effect on the quality of the solution. For all computational variants that we have considered it 
was within the range of several grid cells (typically, eight to ten, see Section 4.4 for further detail). 

Lacunae-based ABCs for (2,2) cell-cente scheme Lacunae-based ABCs for (2,2) node-centered scheme 



Dimensionless time Dimensionless time 


(a) The second-order scheme (2.21) (b) The second-order scheme (2.22) 

Fig. 4.1. Grid convergence study with the lacunae-based ABCs for the wave equation. 


In Figures 4.1 and 4.2 we present the results 
Lacunae-based ABCs for (2,4) node-centered scheme of the grid convergence study for the wave equa- 
tion integrated with the lacunae-based ABCs over 
a long time interval. The errors are evaluated on 
the interior domain S(t) in the maximum norm. 
Computational setup completely follows that of 
Section 2.5, where we applied the lacunae-based 
algorithm to the original problem directly, see Fig- 
ures 2.2 and 2.3. Namely, each scheme was im- 
plemented on a sequence of three grids: 64 x 128, 
128 x 256, and 256 x 512; all geometric parameters, 
grid sizes, parameters that control the partition 
(2.7), as well as the actual exact solution (2.19), 
(2.20), against which we compare our numerical 
Fig. 4.2. Same as Fig. 4-1 for the fourth-order scheme (2.23). results, were taken exactly the same as before. We 

only note that in this case partition (2.7) applies to 
the artificial near-boundary sources needed for constructing of the ABCs, rather than the original right-hand 
side / that drives equation (2.15). 

An obvious observation which is easy to make is that Figures 4.1(a), 4.1(b), and 4.2, look practically 
indistinguishable from Figures 2.2(a), 2.2(b), and 2.3, respectively. In other words, the actual levels of the 
error on corresponding grids are essentially the same. As such, we conclude that in this most simple case the 
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introduction of the ABCs makes the outer boundary of the computational domain completely transparent 
for all the outgoing waves. This is equivalent to saying that the external boundary generates no reflection 
or alternatively, that any imperfections associated with the treatment of the outer boundary can always be 
kept on or below the level of the truncation error pertinent to the interior discretization. In this sense, the 
discrete lacunae-based ABCs that we have constructed can be regarded as an ideal closure of the interior 
finite-difference scheme. Experimentally, this is corroborated by the fact of non-deteriorating convergence of 
the scheme with the theoretically prescribed rate to the specially constructed exact solution of wave-radiation 
type on the computational domain S(t). 

4.2. Nonuniform Motion of the Source. The exact solution (2.19) of the wave equation driven by 
a point source in straightforward uniform motion was obtained in Section 2.5 using Lorentz’ transform. The 
same solution could, of course, have been obtained by directly applying the Kirchhoff integral (2.3). The 
integration would require calculating explicitly the location, at which the trajectory of the source intersects 
the characteristic cone, and would lead to the same analytic result through a somewhat more complex 
derivation. In this section, we consider a somewhat more complex case of straightforward nonuniform (i.e., 
accelerated) motion of the source. 11 The Lorentz transform will obviously not apply in this case, but the 
Kirchhoff integral can still be used for obtaining the exact solution. However, for the general accelerated 
motion it may be impossible to analytically find the intersection of the characteristic cone with the source 
trajectory. Basically, it will require numerical computation, thus making the resulting exact solution only 
“semi-analytic.” As such, to analyze the case of the accelerated motion, we have rather chosen a full- 
fledged numerical approach. First we calculate a fine-grid reference solution using the original lacunae-based 
algorithm of Section 2, and then compare against it the solutions obtained with the ABCs on coarser grids. 

Taking the value of the speed of sound to be equal to one, c — 1, we consider the following law of 
accelerated motion for the center of the domain S(t): 

r — 0, z = Zo(£) — kt + k(cost — 1), 

where k — 0.1. Keeping the values of all geometric parameters (sizes of the domains, etc.) the same as 
before, see Section 2.5, we introduce the following excitation for the wave equation (2,15): 

( st ) ' r ( 5 ) ' (‘ + ' p (* - V) • <41) 

where d = 1.8 k = 0.8, f 2 = r 2 + (z — Zo(t)) 2 , and P(-) is the polynomial of 9th degree that decays smoothly 
from 1 to 0 on the interval [0, 1] so that P( m )(0) = P^ 771 ^ 1) = 0 for m = 1,2,3, and 4 (see Section 2.5). The 
function / of (4.1) obviously has four continuous derivatives everywhere with respect to all its arguments. 
Notice, the temporal behavior of this source function has been purposely chosen sufficiently complex; the 
frequency of the magnitude oscillations and that associated with the motion are incommensurable. 

Equation (2.15) driven by the RHS (4.1) was integrated on the fine grid of dimension 512 x 1024 till 
t = 50 ■ djc using the lacunae-based algorithm of Section 2 implemented with the fourth-order scheme (2.23). 
We have chosen here a shorter time interval compared to those w r e have used for previous demonstrations 
(Sections 2.5 and 4.1) so that not to make the computation of the reference solution excessively expensive. 
This interval is still quite sufficient for experimentally judging the convergence, see Figure 4.3 below. Having 

11 In all numerical examples we consider only straightforward motion because its direction has to be aligned with the z 
direction of the cylindrical coordinate system; otherwise, the symmetry will be lost. In general, this is not a limitation. 
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computed the fine-grid reference solution, we then integrated the same problem on a collection of coarser 
grids using the ABCs and compared the results with this reference solution. In so doing, we have employed 
only the two node-centered schemes: The second-order scheme (2.22) and the fourth-order scheme (2.23). 
The reason is that when both the fine-grid solution and the coarse-grid solution are calculated using a node- 
centered scheme, it is very easy to compare them point-wise (e.g., taking every other, every fourth, etc., 
node of the fine grid). In contradistinction to that, if we were to calculate a coarser-grid solution using 
the cell-centered scheme (2.21), then to compare it against the reference solution we would have had to use 
interpolation on the grid. This has a potential of contaminating the results because of the interpolation 
error, therefore we did not perform the aforementioned comparison for scheme (2.21). 


Lacunae-based ABCs for (2,2) node-centered scheme 



Lacunae-based ABCs for (2,4) node-centered scheme 



(a) The second-order scheme (2.22) (b) The fourth-order scheme (2.23) 

Fig. 4.3. Convergence of the solution of equation (2.15), (4-1) obtained with the ABCs to the fine-grid reference solution. 


Lacunae-based (2,2) node-centered scheme 



(a) The second-order scheme (2.22) 

Fig. 4.4. Same as Figure ^.3, but the solution of (2.15), 


Lacunae-based (2,4) node-centered scheme 



Dimensionless lime 


(b) The fourth-order scheme (2.23) 

(4-1) is obtained with the original lacunae-based algorithm. 
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When computing the solution of equation (2.15) driven by non-uniformly moving source (4.1), we, of 
course, employ the computational domain S(t) that traces the motion of the source. As such, the ABCs are 
set on an artificial boundary that performs accelerated motion. On Figure 4.3, we compare the solutions 
obtained with the help of the lacunae-based ABCs with the fine-grid reference solution. Figure 4.3 show’s 
the dependency of the numerical error on dimensionless time for different computational variants. As before, 
w’e use the sequence of three grids: 64 x 128, 128 x 256, and 256 x 512, to demonstrate the convergence. 
Figure 4.3(a) clearly indicates the second-order convergence of scheme (2.22). For scheme (2.23), we observe 
the fourth-order convergence on Figure 4.3(b). 

To assess the performance of the boundary conditions, we have computed the same solution on the same 
collection of coarser grids (64 x 128, 128 x 256, and 256 x 512) with the same schemes (2.22) and (2.23), but 
with no ABCs, rather using the original lacunae-based algorithm of Section 2 (as we did when computing 
the reference solution). On Figure 4.4, w r e compare the results with the exact solution. As expected, scheme 
(2.22) converges uniformly in time with the second order (see Figure 4.4(a)), and scheme (2.23) — with the 
fourth order (see Figure 4.4(b)). 

Comparing Figure 4.3(a) against Figure 4.4(a) we conclude that for the second-order scheme (2.22), 
the introduction of the ABCs again gives rise to no reflection back into the computational domain (i.e., no 
reflection beyond the level of the truncation error in the interior, cf. Section 4.1). As concerns the fourth- 
order scheme (2.23), one can still notice slight differences between the respective curves on Figures 4.3(b) and 
4.4(b). The difference is most visible for the finest grid 256 x 512, less visible for the medium grid 128 x 256, 
and non-existent for the coarsest grid 64 x 128. This indicates that a small amount of reflections due to the 
ABCs may be present in the solution, although the actual elevation of the error on Figure 4.4(b) compared to 
4.4(a) is so low r that we can regard these reflections negligible anyway. Nonetheless, the discrepancy between 
the corresponding curves needs to be accounted for. We attribute it to the higher sensitivity of the fourth- 
order algorithm to the quality of the discrete lacunae. This phenomenon is commented on in Section 4.4; it is 
not of the fundamental nature, the quality of the lacunae can rather be controlled by appropriately choosing 
the parameters of the numerical procedure, more precisely, the shape and smoothness of the multiplier /x. 

Summarizing for the case of accelerated motion, we see that the solution of non-deteriorating quality 
on long time intervals can still be successfully computed using the lacunae-based ABCs. To the best of 
our knowledge, no other ABCs’ methodology available in the literature can handle artificial boundaries of 
domains that move with acceleration while always keeping the reflections on or below T the level of truncation 
error that pertains to a given interior discretization. 

4.3, Variable Speed of Sound. For the last set of numerical experiments that we present in the 
paper, we wanted to pick up a case that would supposedly be more prone to the buildup of numerical error 
inside the computational domain 5(<). At the same time, we wanted to keep the computational setup “in 
line” with the previous experiments, see Sections 2.5, 4.1, and 4.2. In this connection, we notice that all the 
schemes that we have been using for numerical demonstrations previously were of explicit central-difference 
type. Consequently, the associated discretization error w T as of primarily dispersive nature. In the example of 
the current section, we wall artificially increase the numerical dispersion inside the computational domain S(t) 
and experimentally assess the performance of the combined methodology (interior scheme and the ABCs) for 
this case. We should note, how r ever, that the use of the ABCs is certainly not limited to the aforementioned 
class of the schemes. 

It is known that numerical dispersion for central-difference schemes is more visible for more “suboptimal” 
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Courant numbers. In other words, the further below the stability limit the Courant number is, the more 
dispersive the numerical waves become. In particular, it is easy to see that the one-dimensional second- 
order scheme (3.2b) is exact, and simply reduces to pure propagation along the characteristics, when the 
Courant number ^ is equal to 1. Reducing this number will introduce dispersion of numerical waves. (Of 
course, the convergence of the scheme still implies that the phase shift for every given frequency will become 
smaller as the grid sizes becomes smaller.) The analysis of the one-dimensional case also indicates that in 
multi-dimensional settings numerical dispersion is unavoidable. This is easy to understand already from 
the following qualitative consideration: To guarantee stability for all the waves propagating at an angle 
with respect to the grid lines one has to choose a smaller Courant number, which will necessarily appear 
subopt imal for those waves that propagate along the grid lines. 

As has been mentioned, our intention now is to increase the numerical dispersion inside the computational 
domain and subsequently test the performance of the combined algorithm. To do that, we continuously 
reduce the speed of sound c from the peripheral areas of S(t) to its center. As stability across the entire 
domain will still be limited by the maximum speed of sound, the corresponding Courant number near the 
center will be suboptimal. This will imply higher levels of dispersion closer to the domain center. This will 
also mean that any perturbation that originates in S(t) will stay inside the domain longer compared to the 
previously analyzed cases of constant c. The explanation is obvious — the interior speed of propagation is 
lower. Consequently, we may expect that every particular wave will accumulate more error before it leaves 
the domain S{t). 

We emphasize that accurate quantification of the aforementioned phenomena is not of the central interest 
for discussion in the current paper. Therefore, we do not attempt to quantify the foregoing considerations, 
especially as it may appear technically difficult in any non- trivial setting. However, even on the level of qual- 
itative understanding of the mechanisms of numerical dispersion, it is certainly of interest to experimentally 
assess the performance of the scheme with the lacunae-based ABCs for the case of variable c. 

For our actual computations, we have chosen the following law of variation for the speed of sound inside 
the computational domain S(t): 


c 2 = c 2 


1-Po-P 


( 5 ))- 


(4.2) 


where c is the original constant speed of sound in the far field, d denotes the diameter of the computational 
domain as before, and the polynomial P(-), as well as the quantities f and have been introduced in 
Section 2.5. The constant P 0 in expression (4.2) determines the extent of reduction in the speed of sound at 
the center of S(t ); and we have tried two specific values: Pq = 0.9 and P 0 == 0.99, in our simulations. The 
solution that we have been computing in the case of variable speed of sound is the same traveling-source 
solution (2.19b), (2.20) that we analyzed before. However, instead of equation (2.15) we are now solving 
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where c 2 is defined by (4.2). Substituting ip(r, z , t) of (2.20) into the left hand side of equation (4.3) we obtain 
the source term /(r, z,t) that will obviously differ from that used previously in Sections 2.5 and 4.1. This 
new source term is still compactly supported on the domain S(t) C K 3 for all times, and it now drives the 
solution (2.20) on the entire space. As concerns the methodology for setting the ABCs, it remains exactly 
the same as before. Indeed, we point out that the variation of the speed of sound pertains to the original 
interior problem only. And the auxiliary problem that we solve for the purpose of setting the ABCs is by 
definition formulated with the constant speed of sound throughout its entire domain. 
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Lacunae-based ABCs for (2,2) cell-centered scheme 

Variable speed of sound inside the computational domain 
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Lacunae-based ABCs for (2,2) node-centered scheme 

Variable speed of sound inside the computational domain 
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(a) The second-order scheme (2.21) 


(b) The second-order scheme (2.22) 


Fig. 4.5. Convergence of the solution of equation (4-3), (4-2) w *tk the ABCs to the exact solution (2.20) for Pq = 0.9. 


Lacunae-based ABCs for (2,4) node-centered scheme 
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In Figures 4.5 and 4.6 we present the results 
of the convergence study for the case Pq = 0.9 
(see formula (4.2)) on the same sequence of three 
grids that we have used previously: 64 x 128, 
128 x 256, and 256 x 512. From Figures 4.5(a) 
and (b) we conclude, as before, that the algorithm 
converges with the second order for both schemes 
(2.21) and (2.22); and on Figure 4.6 we again 
observe the fourth-order convergence of scheme 
(2.23). The convergence obviously does not dete- 
riorate as the time elapses, at least till dimension- 
less time reaches the moment 200 d/c when we 
stop the computation. This shows that similarly 
to the previous cases of constant c, the proposed 
numerical procedure in the case of variable speed 
of sound is still capable of providing the solution of 
non-decreasing quality. Note, however, that as the original lacunae-based algorithm cannot be applied to the 
equation with variable c, we cannot directly compare in this case numerical results obtained with the ABCs 
against those obtained with no ABCs as we did before (see, e.g., comparison of the results on Figure 4.3 
with those on Figure 4.4 in Section 4.2). As such, in making a conclusion that the ABCs in this case perform 
practically as well as they did in the previous cases, we rely on the foregoing experimental observation of 
temporally uniform convergence, as well as on the fact that the actual error levels on Figures 4.5 and 4.6 are 
only slightly higher than the respective levels on Figures 4.1 and 4.2. This is expected, because the results 
on Figures 4.1 and 4.2 correspond to numerically reproducing the exact same solution <p(r, z,t) of (2.20) 
using the ABCs but applying them to the original constant-coefficient wave equation in the interior. 

Similar conclusions as to the convergence and quality of the numerical solution can be drawn for the case 


Fig. 4.6. Same as Fig. 4-5 for the fourth- order scheme (2.23). 
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P 0 = 0.99 from looking at Figures 4.7(a) and (b) and Figure 4.8. From the qualitative considerations above 
it follows that this case is supposed to be “tougher” to compute, because the reduction in the speed of sound 
is more significant. In practice, this is manifested by noticeably more oscillatory error profiles, although we 
still clearly see that there is no deterioration of the solution in the long run. Besides, the actual levels of the 
error are somewhat higher compared to the corresponding curves on Figures 4.5(a) and (b) and Figure 4.6. 
This is also expected because the numerical dispersion inside S(t) is supposed to be higher as well. 


Lacunae-based ABCs for (2,2) cell-centered scheme Lacunae-based ABCs for (2,2) node-centered scheme 
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(a) The second-order scheme (2.21) 


( b ) The second-order scheme (2.22) 


Fig. 4.7. Convergence of the solution of equation (4-3), (4-2) to the exact solution (2.20) for Pq = 0.99. 


Lacunae-based ABCs for (2,4) node-centered scheme 
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4.4. Implementation Notes. The forego- 
ing algorithm of lacunae-based ABCs has several 
parameters that need to be tuned appropriately so 
that to obtain the best possible results. Most of 
the flexibility associated with the algorithm resides 
in constructing the multipliers and artificial near- 
boundary sources needed for computing the ABCs 
(Section 3), as well as in choosing the parameters 
of the lacunae-based integration (Section 2). We 
have not yet conducted a comprehensive study of 
how the corresponding parameters affect the nu- 
merical procedure and as such, will only outline 
here some general trends. 


As has been mentioned before, the multiplier 
Fig. 4.8. Same as Fig. 4.7 for the fourth-order scheme (2.23). has to be smooth in the transition region, see Fig- 
ures 3.2 and 3.3. Otherwise, lacunae of the contin- 
uous solution will not be reproduced sufficiently accurately in the discrete solution of the auxiliary problem 
(essentially because the scheme will lose consistency, see discussion in the end of Section 3.2). In most 
of our computations we have used an algebraic polynomial function with four continuous derivatives for 
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multiplier, and the extent of the transition region was about ten grid cells. This has always been sufficient 
for the second-order schemes. In other words, we could always obtain the temporally uniform second-order 
convergence with these settings, although we did not check, for example, whether or not it w T as possible 
to further reduce the extent of the transition region. As concerns the fourth-order convergence, we did see 
situations when the previous settings turned out somewhat insufficient (e.g., in Section 4.2). This occurred 
mostly w T hen going from the medium grid 128 x 256 to the finest grid 256 x 512. To maintain the convergence 
rate in this case, w T e had to use a wider transition region (fifteen cells) and/or a smoother multiplier (five 
continuous derivatives). This indicates that in general the algorithm of lacunae-based ABCs is sensitive to 
the smoothness of the multiplier, as it is supposed to be. However, this sensitivity does not actually manifest 
itself before the error reaches sufficiently low levels. As such, in practical computing one will most likely be 
able to use rather narrow 5 * 7 transition regions, as well as multipliers with limited smoothness. 

As concerns choosing the parameters of the lacunae-based integration (see Section 2), there is at least one 
important observation that has been made experimentally. In theory, the contribution of a given fragment 
of the RHS can be subtracted from the overall solution as soon as the time interval T; nt has elapsed since 
its inception. In practice, it has been found useful to introduce the so-called aft front time gap 5, i.e., allow 7 
for a little extra time for the w r aves to propagate outw r ard. This implies choosing a somewffiat larger value 
for the period Z compared to the necessary minimum given by (2.6) so that by the time of subtraction, 
w T hich is Ti n t + <5 > Ti nt , the reflected w T aves will not have started re-entering the domain S(t) yet. Moreover, 
we can choose to introduce the actual front time gap as w r ell, i.e., increase Z even further, so that by the 
time of subtraction of a given contribution the corresponding reflected waves still be at a (small) distance 
from S(t) rather than right next to its boundary. Experimentally, we have found that the aft front time gap 
affects the quality of the solution stronger than the actual front time gap. However, the quantity 5 in all 
our simulations w r as sufficiently small anyw r ay, about 4% of Ti nt , and most likely it could be reduced even 
further. 

5. Conclusions. We have constructed and tested the algorithm for setting highly-accurate global arti- 

ficial boundary conditions in the problems of time-dependent wave propagation. The key building block of 
the new ABCs is a special non-deteriorating numerical procedure that has been developed previously for the 
long-term integration of wave- radiation problems. The latter procedure is based on the presence of lacunae 

(aft fronts of the waves) in the three-dimensional wave- type solutions. The resulting lacunae-based ABCs 
are obtained directly for the discrete formulation of the problem and can complement any consistent and 
stable finite-difference scheme. In so doing, neither a rational approximation of non-reflecting kernels, nor 
discretization of the continuous boundary conditions is required. The extent of temporal nonlocality of the 
new ABCs appears fixed and limited, and this is not a result of any approximation but rather a direct con- 
sequence of the fundamental properties of the solution. The proposed ABCs can handle artificial boundaries 
of irregular shape on regular grids with no fitting/ adaptation needed. Besides, they possess a unique capa- 
bility of being able to handle boundaries of moving computational domains, including the case of accelerated 
motion. We have conducted a series of numerical experiments that would corroborate the theoretical design 
properties of the algorithm. The experiments included computation of unsteady wave-radiation solutions 
over long time intervals. In all our experiments the ABCs could always keep the level of reflections from the 
artificial boundary on or below the level of truncation error for the interior discretization for as long as the 
computation w^as run. Besides the classical wave equation that we have analyzed in the current paper, the 
proposed technique may find applications in computational acoustics and computational electromagnetics. 
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